Collective motion of self-propelled particles: kinetic phase transition in one dimension 
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We demonstrate that a system of self-propelled particles (SPP) exhibits spontaneous symmetry 
breaking and self-organization in one dimension, in contrast with previous analytical predictions. To 
explain this surprising result we derive a new continuum theory that can account for the development 
of the symmetry broken state and belongs to the same universality class as the discrete SPP model. 



The transport properties of systems consisting of self- 
propelled particles (SPP) have generated much attention 
lately ||l|-0|. This interest has been largely motivated 
by analogous processes taking place in numerous biolog- 
ical phenomena (e.g., bacterial migration on surfaces pi, 
flocking of birds, fish, quadrupeds ||g] , correlated motion 
of ants Q and pedestrians ijl^l ) as well as in various other 
systems, including driven granular materials |ll| , |l2| and 
traffic models [n3[ . The models describing these phenom- 
ena are distinctively non-equilibrium, exhibiting kinetic 
phase transitions and self-organization, of particular in- 
terest from the point of view of modern statistical me- 
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In the simplest version of the SPP model M - intro- 
duced to study collective biological motion - each par- 
ticle's velocity is set to a fixed magnitude, vq. The in- 
teraction with the neighboring particles changes only the 
direction of motion: the particles tend to align their ori- 
entation to the local average velocity. Numerical simu- 
lations in 2D provided evidence of a second order phase 
transition p5[ between an ordered phase in which the 
mean velocity of the entire system, (u), is nonzero and 
a disordered phase with {v) = 0, as the strength of the 
noise is increased or the density of the particles decreased. 

This SPP model is similar to the XY model of classi- 
cal magnetic spins because the velocity of the particles, 
like the local spin of the XY model, has fixed length and 
continuous rotational symmetry. In the vq — Q and low 
noise limit the model reduces exactly to a Monte-Carlo 
dynamics of the XY model. Since the XY model does 
not exhibit a long-range ordered phase at temperatures 
T > [|l^, the ordered state observed in ||l| is surpris- 
ing. To explain this discrepancy. Toner and Tu (TT) 
PI proposed a continuum theory that included in a self- 
consistent way the non-equilibrium effects as well. They 
have shown that their model is different from the XY 
model for d < A and found an ordered phase in d = 2 
[p7|. While TT could account for the first time for the 
ordered phase in 2D SPP models, their theory does not 



allow for an ordered phase for d = 1. 

In this paper we demonstrate that kinetic phase tran- 
sition and ordering takes place in ID as well. This result, 
not foreseen by the existing analytical approaches, moti- 
vated us to introduce a new continuum theory describing 
the SPP model in arbitrary dimensions. Linear stability 
analysis shows that indeed in ID the continuum model 
exhibits ordering. Furthermore, numerical investigations 
indicate that the continuum theory and the discrete SPP 
model belong to the same universality class. 

The ID SPP model — Let us consider N off-lattice 
particles along a line of length L. The particles are char- 
acterized by their coordinate Xi and dimensionless veloc- 
ity Ui updated as 



Xi{t + 1) = Xi(t) + voUi{t), 
u^{t + l)^G({u{t)),) +e:. 



(1) 



The local average velocity {u)i for the ith particle is cal- 
culated over the particles located in the interval [xi — 
A,Xi + A], where we fix A = 1. The function G incorpo- 
rates both the propulsion and friction forces which set the 
velocity in average to a prescribed value vq: G{u) > u 
for w < 1 and G{u) < u for u > 1 iQ. The distribu- 
tion function P(x = f ) of the noise ^,; is uniform in the 
interval [— r;/2, 77/2]. 

Keeping uq constant (uq = 0.1), the adjustable con- 
trol parameters of the model are the average density of 
the particles, p — N/L, and the noise ainplitude r]. We 
implemented one of the simplest choices p9| for G, 



G{u) 



{u + l)/2 forw>0 
(w-l)/2 forw<0. 



(2) 



and applied random initial and periodic boundary con- 
ditions. 

In Fig. 1 we show the time evolution of the model for 
■q = 2.0. In a short time the system reaches an ordered 
state, characterized by a spontaneous broken symmetry 
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and clustering of the particles. In contrast, for larger 
values of ry a disordered velocity field can be found. 

Scaling and exponents — To capture quantitatively the 
transition from an ordered to a disordered state, in Fig. 
2a we plot the order parameter cf) = (u) vs rj for various 
p. As in the two dimensions 0;^ql, the ordered phase 
emerges through a second order phase transition. Near 
the critical noise amplitude, ric{p, L), which separates the 
ordered from the disordered phase, (j) vanishes as 



(t>{v,p) 




for rj < i]c{p,L) 
for rj > i]c (/O, L) 



(3) 



which finding is supported by the increasing scaling 
regime with the system size L (Fig. ^b) and by the con- 
vergence of rjcip, L) to a non-zero rjc{p, oo) value. We 
find that /3 — 0.60 ± 0.05, which is different from both 
the mean-field value 1/2 |^ and /3 = 0.42 ± 0.03 found 
ind=2 |l|. 

Fig. 2a also shows that the various 0(?7, p) curves 
can be collapsed onto a single function (J)q{x), where 
X = r]/r]c{p), just like in d = 2. As shown in Q, the 
consequence of this fact is that near the critical density 
the order parameter vanishes as 



Hv,p) 




ior p > pc{ri, L) ^ ^^-^ 
for p < Pc{v,L) 



with /3' = p. These results can be summarized in the 
p — rj phase diagram shown in Fig. 2c. We also find that 
the critical line, T]c{p), follows 



Vcip) 



(5) 



with K = 0.25 ±0.05. 

While the above numerical results demonstrate the ex- 
istence of the phase transition in one dimension and pro- 
vide numerical values for the scaling exponents /?, (3' and 
K, the origin of these values is unclear at this point. In 
particular, the emergence of the ordered phase in ID is 
not predicted either by the equilibrium theories or by 
the TT model. To overcome this limitation next we 
introduce and investigate a set of continuum equations 
(which can be generalized to any dimension) in terms of 
U{x,t) and g{x,t), where U and g represent the coarse- 
grained dimensionless velocity and density fields, respec- 
tively. Analogous approaches were fruitful in the study 
of a similar SPP system, one-lane traffic flow |M . 

Continuum theory — Let us denote by n(u, x, t)dudx 
the number of particles moving with a velocity in the 
range of [vqu,vq(u + du)] at time t in the [x,x + dx] 
interval. The particle density q{x^ t) is then given as 
g = J ndu, while the local dimensionless average ve- 
locity U{x,t) can be calculated as gll = J nudu. Ac- 
cording to the microscopic rules of the dynamics, in a 
given time interval [t, t -\- t] all particles choose a cer- 
tain velocity v/vq = [G{{u)) + £,] and travel a distance 



VT. Thus, the time development of the ensemble aver- 
age (denoted by overline) of n is governed by the mas- 
ter equation n{u,x,t + t) = g{x',t)p{u\U{x' ,t)), where 
x' — X — vqut and p(u\U) denotes the conditional prob- 
ability of finding a particle with a velocity u when the 
local velocity field U is given. From Eq. (1^) we have 
p{u\U) = P{u — G{{U))). Since n is finite, the actual oc- 
cupation numbers in a given system differ from n. This 
fact can be accounted by adding an intrinsic noise term 
to the master equation as 



n{u, X, t 



q{x, t)p{u\U{x', t)) + u{u, x', t), (6) 



where v has the following properties: (i) V — 0, (ii) due 
to the conservation of the particles J vdu — 0, and (iii) 
since we have a random sampling process, the actual 
values of n satisfy Poisson statistics, i.e., the distribu- 
tion function of v depends on g^ u and U as P(i^) — 
y+^eyij>{~\)/T{u + \ + l), where A = gp{u\U). Thus, 
we have v^ = n. 

Taking the Taylor expansion of n{u, x — vqut, t) up to 
second order in x and integrating Eq. (0) according to 
du, in the vqt <C 1, cr^ = J P{u)u^du 3> 1, £> ^ 1 and 
VQTa^ < 1 limit we obtain 



dtg = -vodxigu) + Ddlg, 



(7) 



where D = v'^ra'^ /2. Note, that the appearance of 
the diffusion term is a consequence of the non-vanishing 
correlation time r. Since ^ p{u\U)udu = G{{U)), in- 
tegrating Eq.(||) according to udu, expanding (U) as 
([/) = U+[dlU + 2[d^U)[d^g)/g\/(!> Q, using (g), we 
arrive at 



dtU = f{U) + fi^dlU + a 



id.U)id,g) 



Q 



c, 



(8) 



where /([/) = {G{U) - U)It, p^ = {dG / dU) / {6t) , 
a — 2p^ and C = Jiyudu/gr. Note, that f{U) is an 
antisymmetric function with f{U) > for < U < 1 
and f{U) < for [/ > 1, C = 0, and p ^ crVgr^. 

At this point we consider Eqs. (Q) and (H) with the 
coefficiens p, a, a, vo and D as the continuum theory 
describing a large class of SPP models. These equations 
differ from both the equilibrium field theories and the 
nonequilibrium system investigated by TT Q . The main 
differences are due to (i) the nonlinear coupling term 
{dxU){dj;g)/ g, and (ii) the statistical properties of the 
noise C. For a = the dynamics of the velocity field U 
is independent of g, and Eq.(^ is equivalent to the $^ 
model describing spin chains, where domains of opposite 
magnetization develop at finite temperatures pO| . 

To study the effect of the nonlinear term in (|8|), we now 
investigate the development of the ordered phase in the 
deterministic case (ct = 0). For a = Eqs. (0) and (||) 
have a set of (meta) stable stationary solutions g* and U* 
describing a "domain wall" separating two regions with 
opposite velocities. Since we can freely translate these 



solutions, we assume t/*(0) — 0. Performing linear sta- 
bility analysis we next show that for certain finite values 
of a the above stationary solutions are unstable. 

Linear Stability Analysis — In the following we make 
use of the fact that the dynamics of g is very slow com- 
pared to that of U. Let us write U in the form of 
U{x,t) = UQ{x,t)+u{x,t), whereUoix,t) = U*{x~^{t)) 
and £,{t) is defined by C/(^(i),t) = 0, i.e. £^{t) defines 
the position of the domain wall. Moreover, we substitute 
f = U — U^ as one of the simplest choices for f{U). Now 
in the u <^ Uq and dxU <^ d^Uo linfit in the moving frame 
x' — X — ^{t) Eq.(^ reads as 

dtu' ^^a+{g- 2)u' + n'^d^u' + ah{x' + S), (9) 

where u'{x') = m(x), [dxg*{x)][dxUa{x)]/ q*{x) == h{x' + 
0, dxU*{x') = a{x') and 5(2;') = {dJ/dU){U*{x')) + 2. 
Fromf7(^(t),t) = Uo{£.{t),t) ^Owegetu'{x' = 0,t) = 0, 
which yields 



Ca(0) = fi^dlu'(Q,t) + ah{(,). 



(10) 



Now (p) and ( [lO[ ) describe the time development of the 
velocity field in terms of u' and ^. 

Fourier transforming (P) and (ll^) we find that the 
short wavelength fluctuations (fc —> cxj) are stabilized by 
the Laplacian term in (H). However, the growth rate A of 
the long wavelength fluctuations, defined by ^ ~ exp(At), 
is determined by the characteristic equation 

(_2 + V-A)(^^-A)-al^^ = 0. (11) 



V2-D/i 



D 



It can be seen that (Oil) has positive solution for large 
enough a. This result means that for a > ac {ac ~ 
2^/2Dii/vo) the domain wall solution U* is unstable, 
hence in this regime the walls disappear and all parti- 
cles move in the same direction, demonstrating that (R) 
and M) predict an ordered phase in one dimension. 

To further confirm the relevance of the continuum the- 
ory to the discrete model (1), we integrate numerically 
(0) and (^). The results of the integration are in excel- 
lent agreement with the results obtained for the discrete 
model and with the linear stability analysis, and can be 
summarized as follows: (i) For the noiseless case (u = 0) 
we find that for a > ac there is an ordered phase, which 
disappear for a < ac, (ii) The ordered phase predicted 
by ( [ill ) is present for the noisy ct > as well; (iii) In- 
creasing a leads to a second order phase transition from 
the ordered to the disordered state. Since a plays the 
role of J] in (|l|) , this transition is equivalent with (3) ob- 
served in the discrete model; (iv) Finally, we measured 
the order parameter as a function of a for various val- 
ues of D. As Fig 2d. illustrates, for small D Eq.(3) with 
/3 = 0.6 provides an excellent fit to the numerical results, 
indicating that the discrete model (1) and the continuum 
theory (M) and (0) belong to the same universality class. 



In conclusion, we showed that the SPP model exhibits 
spontaneous symmetry breaking and ordering in one di- 
mension, a result surprising both in the light of the equi- 
librium spin models, and of the continuum theory inves- 
tigated by TT. We introduced a new continuum theory, 
whose terms are explicitelly derived from the ingredients 
of the discrete model (in contrast with constructing it 
from symmetry arguments). Linear stability analysis in- 
dicates that an ordered phase can develop as the domain 
walls become unstable. While here we limited ourselves 
to the ID case of the continuum theory, (1^) and (pi) can 
be generalized to arbitrary dimensions, thus apply to the 
physically and biologically relevant two and three dimen- 
sional models as well. Since the continuum model inves- 
tigated by us contains the major ingredients of the SPP 
models, such as self-propulsion and local reorientation of 
the velocity, we expect our results to apply to a wide 
variety of systems made of SPP particles ||,||,^,0,|j . 
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FIG. 1. The dynamics of the ID SPP model for L = 300, ry = 2.0 and A'^ = 600. The darker gray level represents higher 
particle density. Note that the particles exhibit clustering and the spontaneous broken symmetry of motion. 



FIG. 2. a: The order parameter (j) vs the noise amplitude normalized by the critical amplitude r]c{p), for L — 1000 and 
various values of p. For -q < rjc{p) the system is in a symmetry-broken state indicated by (/) > 0. b: cf) vanishes as a power-law 
in the vicinity of r]c{p). Note the increasing scaling regime with increasing L. The solid line is a power-law fit with an exponent 
(3 = 0.6, while the dotted line shows the mean-field slope (3 = 1/2 as a comparison, c: Phase diagram in the p — rj plane. The 
critical line follows 77c (p) ~ p". The solid curve represents a fit with k = 1/4. d: The order parameter vs the std deviation of 
the noise normalized by a* — Oc(D = 0), obtained by direct numerical integration of the continuum model for a = 2, p = 1, 
vo = 0.1, p = 1, L = 1000 and various values of D. For D <C 1 0((t) follows a power-law with an exponent /3 — 0.6 (solid line). 
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